Simple Linear Regression Models

STA 101 - Summer I 2022

Raphael Morsomme

Welcome

Announcements

Assignments

Recap of last lecture

Histogram, scatterplot, boxplot
Average, median, variance, sd and IQR; robustness
Frequency, contigency and proportion tables
Barplot, mosaic plot
Effective communication: well-edited figures, $\ge3$ variables (symbols, colors, facets)
R for Data Science - chapters 3 and 7

Outline

simple linear regression model
least-square estimates
model comparison
outliers

Simple linear regression

Regression model

d <- ggplot2::mpg
head(d, n = 4)

# A tibble: 4 x 11
  manufacturer model displ  year   cyl trans      drv     cty   hwy fl    class 
  <chr>        <chr> <dbl> <int> <int> <chr>      <chr> <int> <int> <chr> <chr> 
1 audi         a4      1.8  1999     4 auto(l5)   f        18    29 p     compa~
2 audi         a4      1.8  1999     4 manual(m5) f        21    29 p     compa~
3 audi         a4      2    2008     4 manual(m6) f        20    31 p     compa~
4 audi         a4      2    2008     4 auto(av)   f        21    30 p     compa~

Suppose that the variable hwy (fuel efficiency on the highway) is very expensive to measure.

We decide to estimate it using the other variables. To do so, we will fit a regression model. \[ \text{hwy} \approx \text{model}(\text{other variables}) \]

Simple regression

We expect the variable cty to be a good proxy for hwy.

After all, if a car is efficient in the city, we expect it to also be efficient on the highway! We will therefore consider a simple regression model

\[ \text{hwy} \approx \text{model}(\text{cty}) \]

Simple linear regression

ggplot(d) +
  geom_point(aes(x = cty, y = hwy))

It turns out that the variables cty and hwy are linearly associated.

We therefore decide to use a **simple linear regression* model

\[ \text{hwy} \approx \beta_0 + \beta_1 \text{cty} \]

Parameters

The numbers $\beta_0$ and $\beta_1$ are the two parameters of the model.

We now want to find good values for these parameters.

Let us look at two sets of values:

$\beta_0 = 1$, $\beta_1 = 1.3$ \[ \text{hwy} \approx 1 + 1.3 \text{cty} \]
$\beta_0 = 15$, $\beta_1 = 0.5$ \[ \text{hwy} \approx 15 + 0.5\text{cty} \]

Good model
Bad model

ggplot(d) +
  geom_point(aes(x = cty, y = hwy)) +
  geom_abline(intercept = 1, slope = 1.3)

ggplot(d) +
  geom_point(aes(x = cty, y = hwy)) +
  geom_abline(intercept = 15, slope = 0.5)

Prediction

We can use our models to estimate hwy for new vehicles.

Imagine there is a new vehicle with $\text{cty} = 30$. Instead of measuring its hwy, we decide to use our model to estimate it. Using the “good” model gives the following estimate \begin{align*} & _0 + _1 \ & = 1 + 1.3 * 30\ & = 40\ \end{align}

Is this the true hwy os the new vehicle? No! - this is only an estimate based on the value of the variable cty and our “good” model.

Can we do better? Yes! - Take additional variable into account in the model (e.g. engine size, vehicle age, etc) - Use better values for $\beta_0$ and $\beta_1$.

Group exercise - parameters

What is the prediction for the new vehicle if we use the bad model?
copy andpaste the following piece of code and try different values for the parameters to find a good set of values.

library(tidyverse)
d <- ggplot2::mpg
beta_0 <- 15
beta_1 <- 0.5
ggplot(d) +
  geom_point(aes(x = cty, y = hwy)) +
  geom_abline(intercept = beta_0, slope = beta_1)

05:00

Linear association

A simple linear regression model is only applicable if the relation between the predictor and the response is linear.

If the relation is not linear, the simple linear regression is not suitable.

In this case, we need to model the non-linearity (next lecture).

Group exercise - linear association

Exercise 7.3

03:00

Leat-square estimates

Residuals

Our predictions are only approximate.

Let us represent our prediction with $\widehat{\text{hwy}}$ and the true value with $\text{hw}$
the error we make is $\text{hw} - \widehat{\text{hwy}}$
this is the residual

\[e = \text{hwy} - \widehat{\text{hwy}}\]

Visualizing residuals

Black circles: Observed values (y = hwy)
Pink solid line: Least-squares regression line
Maroon triangles: Predicted values (y = .fitted)
Gray lines: Residuals

Residual plot

d %>%
  mutate(
    hwy_hat = 1 + 1.3 * cty,
    resid   = hwy - hwy_hat
    ) %>%
  ggplot() +
  geom_point(aes(cty, resid)) +
  geom_abline(intercept = 0, slope = 0, col = "red")

Group exercise - residuals

Exercises 7.1, 7.17, 7.19

05:00

Good estimates

We want to choose estimates that give accurate prediction
We want to minimize the residuals

Perhaps the most natural thing to do is to find the values for $\beta_0$ and $\beta_1$ that minimize the sum of absolute residuals \[|e_1|+|e_2|+\dots+|e_n|\] For practical reasons, the sum of squared residuals is more common criterion \[e_1^2+e_2^2+\dots+e_n^2\] ## Why squaring the residuals?

can work by hand (pre-computer)
reflects the assumptions that being off by $4$ is more than twice as bad as being off by $2$
nice mathematical properties
mainstream

Least-square estimates

We find the values for $\beta_0$ and $\beta_1$ that minimize the SSR with the R function lm

lm(hwy ~ cty, data = d)


Call:
lm(formula = hwy ~ cty, data = d)

Coefficients:
(Intercept)          cty  
      0.892        1.337

We denote the least-square estimates with the symbols $\hat{\beta}_0$ and $\hat{\beta}_1$.

In our model, we therefore have \[\hat{\beta}_0 = 0.892, \qquad \hat{\beta}_1 = 1.337\] ##

ggplot(d) +
  geom_point(aes(x = cty, y = hwy)) +
  geom_abline(intercept = 0.892, slope = 1.337)

Reading `R` output

When you fit a model, R can provide the results in different formats.

m <- lm(hwy ~ cty, data = d)

Estimates only
Raw R output
Tidy R output


Call:
lm(formula = hwy ~ cty, data = d)

Coefficients:
(Intercept)          cty  
      0.892        1.337

Printing the model itself provides the least-square estimates. This is sufficient for now.

summary(m)


Call:
lm(formula = hwy ~ cty, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3408 -1.2790  0.0214  1.0338  4.0461 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.89204    0.46895   1.902   0.0584 .  
cty          1.33746    0.02697  49.585   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.752 on 232 degrees of freedom
Multiple R-squared:  0.9138,    Adjusted R-squared:  0.9134 
F-statistic:  2459 on 1 and 232 DF,  p-value: < 2.2e-16

This format contains information that will be useful when we do inference. It is, however, difficult to read.

summary(m) %>%
  tidy() %>%
  kable()

term	estimate	std.error	statistic	p.value
(Intercept)	0.8920411	0.4689457	1.902227	0.05838
cty	1.3374556	0.0269732	49.584698	0.00000

This format contains the information necessary for doing inference and is easy to read. This is the format used the textbook.

Exploring the model

We can access the residuals and prediction (fitted values) with the command augment.

m <- lm(hwy ~ cty, data = d)
augment(m)[,1:4]

# A tibble: 234 x 4
     hwy   cty .fitted .resid
   <int> <int>   <dbl>  <dbl>
 1    29    18    25.0 4.03  
 2    29    21    29.0 0.0214
 3    31    20    27.6 3.36  
 4    30    21    29.0 1.02  
 5    26    16    22.3 3.71  
 6    26    18    25.0 1.03  
 7    27    18    25.0 2.03  
 8    26    18    25.0 1.03  
 9    25    16    22.3 2.71  
10    28    20    27.6 0.359 
# ... with 224 more rows

Group exercise - exploring a model

For the third observation,

What is the value of cty?
What is the predicted value for hwy based on the model?
What is the actual value for hwy?
What is its residual?
Did the model over- or under-predict?

04:00

Intepreting the parameters

the intercept $\hat{\beta}_0$ is the estimate for a car with $\text{cty} = 0$
- meaningless in our case

the slope $\hat{\beta}_1$ is the increase in a our estimate for hwy for each additional unit of
- For each additional unit of cty, we expect hwy to increase by 1.337.

extrapolation
- stick to the range of the data
- the intercept will not always be meaningful

Special case: categorical predictor

Let us create a binary predictor indicating if the car is from $1999$.

d_binary <- mutate(d, year_binary = year == 1999)
head(select(d_binary, hwy, year, year_binary))

# A tibble: 6 x 3
    hwy  year year_binary
  <int> <int> <lgl>      
1    29  1999 TRUE       
2    29  1999 TRUE       
3    31  2008 FALSE      
4    30  2008 FALSE      
5    26  1999 TRUE       
6    26  1999 TRUE

The variable year_binary takes the value $1$ (TRUE) if the car is from $1999$ and $0$ (FALSE) otherwise.

lm(hwy ~ year_binary, data = d_binary)


Call:
lm(formula = hwy ~ year_binary, data = d_binary)

Coefficients:
    (Intercept)  year_binaryTRUE  
       23.45299         -0.02564

The model equation is

\[ \text{hwy} \approx 23.46 - 0.026 \text{year_binary} \]

For a new car from $1999$, the variable year_binary takes the value $1$ and our prediction is

\[ \widehat{\text{hwy}} = 23.46 - 0.026*1 = 23.434 \]

While for a car that not from $1999$, the variable year_binary takes the value $1$ and our prediction is

\[ \widehat{\text{hwy}} = 23.46 - 0.026*0 = 23.46 \]

Model comparison

Alternative model

Let us fit an alternative model using engine size (disp) as a predictor \[ \text{hwy} \approx \beta_0 + \beta_1 \text{displ} \]

ggplot(d) +
  geom_point(aes(displ, hwy))

The least-square estimates for the coefficients are

lm(hwy ~ displ, data = d)


Call:
lm(formula = hwy ~ displ, data = d)

Coefficients:
(Intercept)        displ  
     35.698       -3.531

Note that the slope coefficient is negative; which makes sense since cars with larger engines would tend to be less efficient.

We now have two models. Which is the best?

We could start by looking at the residuals

Comparing residuals

m <- lm(hwy ~ cty, data = d)
m_augment <- augment(m)
ggplot(m_augment) +
  geom_histogram(aes(.resid))

lm(hwy ~ displ, data = d) %>%
  augment %>%
  ggplot() +
  geom_histogram(aes(.resid))

The first model seems to have residuals of smaller magnitude.

But we need a more systematic approach for comparing models:

looking at a plot can be misleading (illusions)
difficult to compare models with similar residuals

SSR

Instead of comparing histograms to identify the model with the smaller residuals, we can compute the SSR (sum of squared residuals!)

small residuals will give a small SSR
large residuals will give a large SSR

\[ SSR = r_1^2+r_2^2+\dots+r^2_n \]

Simply choose the model with the smaller SSR!

📋 The textbook uses the term SSE (sum of squared errors).

$R^2$

The SSR is also useful in describing the strength of the model.

The SST (total sum of squares) is the sum of squared distance to the mean. \[ SST = (x_1 - \bar{x})^2 + (x_2 - \bar{x})^2 + \dots + (x_n - \bar{x})^2 \] It measured the total amount of variability in the data.

Remember the formula for SSR \[ SSR = (x_1 - \hat{x})^2 + (x_2 - \hat{x})^2 + \dots + (x_n - \hat{x})^2 \] is the amount of variability in the data left unexplained by the model.

By analogy, $SST - SSR$ is the amount of variation explained by the model: \[ \text{data} = SST = SSR + (SST-SSR) = \text{residuals} + \text{model} \]

The coefficient of determination $R^2$ measures the proportion of variation in the data that is explained by the model

\[ R^2 = 1-\dfrac{SSR}{SST} = \dfrac{SST-SSR}{SST} \]

Note that $0\le R^2 \le 1$.

A good model will have small residuals, a small SSR and therefore a large $R^2
A great model will have a $R^2$ close to $1$.
A poor model has an $R^2$ close to 0 (SSR close to SST)

Computing $R^2$ in `R`

To compute $R^2$ in R, simply use the command glance.

m_glance <- glance(m)
m_glance

# A tibble: 1 x 12
  r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.914         0.913  1.75     2459. 1.87e-125     1  -462.  931.  941.
# ... with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

m_glance$r.squared

[1] 0.9137752

lm(hwy ~ displ, data = d) %>% glance %>% .[["r.squared"]]

[1] 0.5867867

The model with cty as a predictor has a $R^2$ value of $0.914$, and the model that uses displ has a $R^2$ of $0.59$ (worse).

We would therefore prefer the first model.

Group exercise - Interpretation

Exercise 7.23 (a-d)
Compute $R^2$ “by hand” for the model with parameters $(\beta_0, \beta_1) = (1, 1.3)$. You can use the following R code.

library(tidyverse)
d <- ggplot2::mpg
cty <- d$cty
hwy <- d$hwy

03:00

Outliers

Remember, in a boxplot, outliers are observation far from the bulk of the data

ggplot(d) +
  geom_boxplot(aes(hwy))

In the context of regression model, an outlier is: an observation that falls far from the cloud of points

Group exercise - outlier

Look at the scatterplot for our data. Are there outliers?

ggplot(d) +
  geom_point(aes(cty, hwy))

Exercise 7.25

03:00

Outliers, leverage and influential points

outlier: observation that falls far from the cloud of points
High leverage point: observation that falls horizontally away from the cloud of points
influential point: observation that influences the slope of the regression line

All influential points are high leverage points.
All leverage points are outliers.
(Venn Diagram)

Dealing with outliers

In regression, outliers have the potential to highly influence the least-square estimates - the least-square estimates are not robust to the presence of outliers

Impact of outliers

Let us contaminate the data with an outlier (cty $=10$ and hwy $=1000$) and fit the same regression model.

d_contaminated <- select(d, cty, hwy) %>%
  add_row(cty = 40, hwy = 400) # add outlier
lm(hwy ~ cty, data = d_contaminated)


Call:
lm(formula = hwy ~ cty, data = d_contaminated)

Coefficients:
(Intercept)          cty  
    -26.046        3.013

The value of the slope has tripled!

ggplot(d_contaminated) +
  geom_point(aes(cty, hwy)) +
  geom_abline(intercept = -26.046, slope = 3.013)

The regression line not longer fits the data well.

Recap

simple linear regression model \[ \text{hwy} \approx \beta_0 + \beta_1 \text{cty} \]
residuals
least-square estimates
parameter interpretation
model comparison with $R^2%
outliers

Simple Linear Regression Models

Welcome

Announcements

Recap of last lecture

Outline

Simple linear regression

Regression model

Simple regression

Simple linear regression

Parameters

Prediction

Linear association

Leat-square estimates

Residuals

Visualizing residuals

Residual plot

Good estimates

Least-square estimates

Reading `R` output

Exploring the model

Intepreting the parameters

Special case: categorical predictor

Model comparison

Alternative model

Comparing residuals

SSR

\(R^2\)

Computing \(R^2\) in `R`

Outliers

Outliers, leverage and influential points

Dealing with outliers

Impact of outliers

Recap

Recap

Simple Linear Regression Models

Welcome

Announcements

Recap of last lecture

Outline

Simple linear regression

Regression model

Simple regression

Simple linear regression

Parameters

Prediction

Linear association

Leat-square estimates

Residuals

Visualizing residuals

Residual plot

Good estimates

Least-square estimates

Reading R output

Exploring the model

Intepreting the parameters

Special case: categorical predictor

Model comparison

Alternative model

Comparing residuals

SSR

\(R^2\)

Computing \(R^2\) in R

Outliers

Outliers, leverage and influential points

Dealing with outliers

Impact of outliers

Recap

Recap

Reading `R` output

Computing \(R^2\) in `R`